%设置页面边距（word标准页面）
\documentclass[a4paper]{article}
\usepackage{geometry}
\geometry{a4paper,left=2.7cm,right=2.7cm,top=2.54cm,bottom=2.54cm}
 
%导入ctex包
\usepackage[UTF8,heading=true]{ctex}

\usepackage{adjustbox} % 用于调整表格大小

% 设置段落首行缩进
\setlength{\parindent}{1.5em}
 
%设置摘要格式
\usepackage{abstract}
\setlength{\abstitleskip}{0em}
\setlength{\absleftindent}{0pt}
\setlength{\absrightindent}{0pt}
\setlength{\absparsep}{0em}
\renewcommand{\abstractname}{\textbf{\zihao{4}{摘要}}}
\renewcommand{\abstracttextfont}{\zihao{-4}} %设置摘要正文字号
 
%设置页眉和页脚，只显示页脚居中页码
\usepackage{fancyhdr}
\pagestyle{plain}
 
%调用数学公式包
\usepackage{amssymb}
\usepackage{amsmath}
 
%调用浮动包
\usepackage{float}
\usepackage{subfig}
\captionsetup[figure]{labelsep=space} %去除图标题的冒号
\captionsetup[table]{labelsep=space} %去除表格标题的冒号
%设置标题格式
\ctexset {
	%设置一级标题的格式
	section = {
		name={,、},
		number=\chinese{section}, %设置中文版的标题
		aftername=,
	},
	%设置三级标题的格式
	subsubsection = {
		format += \zihao{-4} % 设置三级标题的字号
	}
}
 
 
%使得英文字体都为Time NewTown
%\usepackage{times}
 
%图片包导入
\usepackage{graphicx}
\graphicspath{{image/}} %图片在当前目录下的image目录
 
%参考文献包引入
\usepackage{cite}
\usepackage[numbers,sort&compress]{natbib}
 
%代码格式
\usepackage{listings}
\usepackage{graphicx}%写入python代码
\usepackage{pythonhighlight}%python代码高亮显示
\lstset{
	%numbers=left, %设置行号位置
	%	numberstyle=\tiny, %设置行号大小
	keywordstyle=\color{red}, %设置关键字颜色
	commentstyle=\color[cmyk]{1,0,1,0}, %设置注释颜色
	escapeinside=``, %逃逸字符(1左面的键)，用于显示中文
	breaklines, %自动折行
	extendedchars=false, %解决代码跨页时，章节标题，页眉等汉字不显示的问题
	xleftmargin=1em,xrightmargin=1em, aboveskip=1em, %设置边距
	tabsize=4, %设置tab空格数
	showspaces=false %不显示空格
}


 
 
\renewcommand{\refname}{}
 
%item包
\usepackage{enumitem}
 
%表格加粗
\usepackage{booktabs}
 
%设置表格间距
\usepackage{caption}
 
%允许表格长跨页
\usepackage{longtable}
 
%伪代码用到的宏包
\usepackage{algorithmic}
\usepackage{algorithm}
 
%正文区
\title{农作物的种植策略} 
\date{} %不显示日期
 
%文档
\begin{document}
	\maketitle
	\vspace{-6em} %设置摘要与标题的间距
	\zihao{-4} %设置正文字号
	%摘要部分
	\begin{abstract}
			随着全国耕地保护政策的不断执行，中国可供耕地面积连续三年不断增长，但是人均耕地面积仍然不足。在这种背
		景下，寻找最优粮食种植方案，以求在有限的耕地资源中实现收益最大化变得尤为重要。通过科学合理的种植规划
		和技术改进，不仅可以提高粮食产量，还能确保耕地资源的可持续利用，进而保障国家粮食安全。本文针对华北某
		山村，农作物种类和耕作类型两个角度对该问题进行了研究 \par
			针对问题一，本文首先站在2023年的角度上模拟决策过程，从耕地面积、重茬种植、三年种豆和水浇地四大方向设
		置约束条件，建立了以最优耕地方式为目标函数的单目标优化模型。\par
			针对问题二，考虑到农作物的预期销售量、亩产量、种植成本和销售价格存在波动，因此考虑选择利用蒙特卡洛模
		拟的方式解决，将最大化期望收益与收益的方差相结合，最终建立目标函数。 \par
			针对问题三，以上一问的结果为基础，综合考虑现实条件中菜品收益存在替代性和互补性等因素，使用协方差矩阵
		描述预期销售量、销售价格与耕作成本之间的相性。在结合历史数据进行回溯分析，最终得出目标方案。 \\
		%关键词（上文最后一段要用“\\”换行）
		\newline
		\noindent{\textbf{关键词：} \textbf{单目标优化模型}\quad   \textbf{蒙特卡洛模拟}\quad \textbf{替代性和互补性} \quad \textbf{相关性}\quad \textbf{优化模型}} 
	\end{abstract}
	
	\clearpage %换页
	
	%正文部分
	%Part one
	\section{问题背景与重述}
	\subsection{问题背景}
		自党的二十大以来，我国在土地资源的节约和集约利用方面取得了显著的成效。最新的统计数据显示，全国耕地面积
	已经连续三年呈现增长趋势，目前全国耕地总面积达到了19.18亿亩。尽管这一数字看似庞大，但由于我国人口众多，
	人均耕地面积仍然相对较小，仅为1.59亩。这一数据表明，虽然总体耕地面积有所增加，但是人均耕地面积稀缺的问
	题依然存在，这给我国农业的可持续发展带来了挑战。\par
		在此背景下，计划耕作的重要性愈发凸显。通过实施科学合理的耕作计划，不仅可以提高土地的利用率，还能有效缓
	解因人均耕地面积不足所带来的压力。计划耕作要求根据各地不同的自然条件、土壤特性及市场需求，制定出符合实
	际的种植方案，以达到最佳的经济效益和生态效益。\par
		优化作物的种植结构布局不仅是现代农业发展的必然要求，也是实现农业可持续发展的关键环节。通过不断改进种植
	技术和管理方法，不仅可以提升农业生产能力，还能促进人与自然和谐共生，为未来的粮食安全和社会经济稳定打下坚实
	基础。
		
	
	\subsection{问题重述}
	\begin{enumerate}[itemindent=0.5cm ]
		\item 问题一 \par
		\hspace{2em}问题一可总结为 “农作物种植方案优化”。假定各类农作物未来预期销售量、种植成本等相对2023年保持
		稳定，需要在未来种植过程中，合理利用乡村现有耕地和大棚资源，分别在超产滞销和超产降价销售的情况下，达
		到最大利益。
		\item 问题二 \par
		\hspace{2em}问题二可总结为“综合不确定性与风险的种植策略优化”。在制定2024年至2030年的农作物种植方案时，必
		须充分考虑到各种不确定性因素对乡村农业生产的潜在影响。首先，农作物的预期销售量存在波动，如小麦和玉
		米的销售量预计每年增长5\%至10\%，而其他作物的销量可能有±5\%的变化。此外，亩产量由于气候等自然条件
		的影响，每年可能会有±10\%的波动，进一步增加了种植规划的复杂性。与此同时，种植成本预计每年增长5\%，
		这将直接影响到生产利润。对于销售价格，粮食类作物价格相对稳定，而蔬菜价格每年上涨约5\%，食用菌类价
		格则可能逐年下降。基于这些不确定性因素以及可能的市场和气候风险，需要通过合理的数学模型，优化作物
		种植组合，确保在波动条件下的收益最大化，并规避潜在的种植风险。
		\item 问题三 \par
		\hspace{2em}问题三可总结为“考虑作物相关性的最优种植策略优化”，在问题2的基础上，进一步考虑农作物之间的可
		替代性与互补性，以及预期销售量、销售价格、种植成本之间可能存在的相关性，制定2024年至2030年乡村农
		作物的最优种植策略，并与问题2中的结果利用蒙特卡洛模拟方法进行对比分析。\par

	\end{enumerate}
	
	%Part Two
	\section{问题分析}
	\subsection{问题一的分析}
	问题一要求为乡村推荐2024$\sim$2030年7年范围内分别在农作物滞销或将过多部分按照2023年度售价的50\%出售两
	种情况下寻找最优种植方案。由于题目中并未给出以往任何年度的市场蔬菜容纳量，不妨在此假设2023年度所有蔬菜
	全部刚好售出，以寻找问题所求解。考虑到实际问题，由于耕地种植面积不宜太小，将模型简化为可用耕地最小为1
	亩，若当前剩余耕地面积小于1亩，则不进行考虑。
	\subsection{问题二的分析}
	问题二要求在全面考量各种农作物的预期销量、单位面积产量、种植成本及销售价格的不确定性，并考虑到潜在的种植风
	险的情况下，为该乡村制定出从2024年至2030年的最佳种植方案。智慧大棚部分蔬菜在统一区域内会种植两次，按照先前
	规划，相当于将对应面积扩大至2倍，以减少运算成本，但由于计算机代码处理逻辑过于复杂，与原有规划相违背，考虑到
	附件2中所述：智慧大棚第一季可种植的蔬菜作物及其亩产量、种植成本和销售价格均与普通大棚相同且这部分蔬菜价值
	在整体蔬菜中所占比例不高，后期人工根据普通大棚的种植方案对数据进行拟合处理。


	\subsection{问题三的分析}
	
	%Part Three
	\section{模型假设}
	%假设的列表
	\begin{enumerate} 
		\item 假设2023年度所产蔬菜恰好全部售出，即以2023年度蔬菜产出量为市场容纳量。
		\item 假设未来2024到2030年间种植蔬菜不受自然灾害的影响。
		\item 假设每种作物的种植成本无波动。
		\item 假设所有农作物的销售价格与市场需求紧密相关，且市场需求在2024到2030年间保持稳定或按照已知的趋势变化。
		\item 假设小麦和玉米每年增长率未为5\%至10\%, 其他作物预期销售量变化为 $\pm 5\%$。
		\item 假设耕地亩产量每年变化为 $\pm 10\%$。
	\end{enumerate}
	
	%Part Four
	\section{符号说明}
	%浮动体表格，使用table实现
	\begin{table}[H] %[h]表示在此处添加浮动体，默认为tbf，即页面顶部、底部和空白处添加
		\captionsetup{skip=4pt} % 设置标题与表格的间距为4pt
		\centering
		\setlength{\arrayrulewidth}{2pt} % 设置表格线条宽度为1pt
		\begin{tabular}{cc} %c表示居中，l表示左对齐，r表示右对齐，中间添加“|”表示竖线
			\hline
			\makebox[0.15\textwidth][c]{符号} & \makebox[0.6\textwidth][c]{说明}  \\ 
			\hline\\
			$x_{ijk}$ & 表示第 i 块地中 j 号作物在当前的种植面积， k 为植物种植周期标记 $k\in \textbf{(1, 2)}$ \\\\
			$x_{ijk}^{t}$ & 表示第 i 块地中 j 号作物在第 t 年的种植面积, k 为植物种植周期标记 $k\in \textbf{(1, 2)}$  \\\\
			$P_{ij}$ & 表示第 i 块地中是否可种植 j 号作物 $P\in (0, 1)$ \\\\
			$C_{ij}$ & 表示第 i 块作地中 j 号作物的成本（元/亩）  \\\\
			$Q_{ij}$ & 表示第 i 块地中 j 号地的产量  \\\\
			$Z_{i}$ & 表示第 i 块地的面积 \\\\
			$pri_{ij}$ & 表示第 i 块地中 j 号作物的售价（元）   \\\\
			$need_{j}$ & 表示第 j 号作物在2023年度需求量 \\\\
			\hline
		\end{tabular}
		% \hline是横线，采用\makebox设置列宽
	\end{table}
		
	%Part Five
	\section{模型的建立与求解}
	\subsection{问题一模型的建立与求解}
	本题要求在各农作物预期销售量、种植成本等参数相对2023年保持问稳定的前提下，为该乡村寻找2024 $\sim$ 2030年间农作物的最优种植方案。
	\subsubsection{模型的建立}
	\textbf{1. 目标函数的确定}\par
	通过对问题一的分析，该问题属于单目标优化问题。“利润 = 所有货物的售价之和 - 成本 - 滞销造成的浪费
	（或 + 滞销部分售价的50\%）”，由此建立目标函数考虑可能干扰因素和造成影响逻辑后，便可建立优化
	模型，对乡村耕地分配进行决策的过程进行描述。\par
	分别采用变量 i 、j 代表耕地编号和所种植作物编号，用$x$代表作物种植面积、Q代表产量、$pri$代表对应作物单价、$C$代表对应作物每亩所
	需成本、$need$代表对应作物市场需求量。如此规定后，可确定目标为单目标优化函数， 并将最大利润$\max$表示为：\par
	情况一：
	\begin{equation}
		\max \; Y=\sum_{i=1}^{54}\sum_{j = 1}^{41}(x_{ijk}^{t} \cdot Q_{ij} \cdot pri_{ij}-x_{ij}^{t} \cdot C_{ij})-\sum_{j=1}^{41}max\left \{  \sum_{i=1}^{54}{x_{ij}^t \cdot Q_{ij}}-need_j,0\right \} \cdot pri_{ij}
		\label{formula_1} %公式一引用标签， 通过 \eqref{formula_1}引用
	\end{equation}

	\begin{samepage}
	情况二：
	\begin{equation}
		\max \; Y=\sum_{i=1}^{54}\sum_{j = 1}^{41}(x_{ijk}^{t} \cdot Q_{ij} \cdot pri_{ij}-x_{ij}^{t} \cdot C_{ij})-\sum_{j=1}^{41}max\left \{  \sum_{i=1}^{54}{x_{ij}^{t} \cdot Q_{ij}}-need_j,0\right \} \cdot pri_{ij} \cdot 50\%
		\label{formula_2} %公式一引用标签， 通过 \eqref{formula_1}引用
	\end{equation}
	\end{samepage}

	\begin{samepage}
		\textbf{2. 约束条件的确定}\par
			\begin{itemize}
				\item[(a)] 单位土地面积约束\par
				在目标函数$\eqref{formula_1}$中给出了$x_{ij}$变量，但未给出具体约束形式。故以下内容将具体围绕这一重要参数展开。
				由实际情况可知，无论当前作物所需耕地面积多大，都不可能超过每块土地所能支持面积的最大值，可确定单位土地面积约束，即以下约束：
					\begin{equation}
						\begin{split}
							\sum_{j = 1}^{41} x_{ijk}^t\leq Z_{i}, 
							 & 
							k \in \{1, 2\}
						\end{split}
						\label{formula_3} %公式二引用标签， 通过 \eqref{formula_2}引用
					\end{equation}
				\item[(b)] 重茬种植约束\par
				再有，重茬种植时，在当前耕地允许种植该作物的条件下，所占用土地面积也不能超过当前剩余的最大土地面积，根据此条件，可确立重茬种植约束如下：
					\begin{equation}
						\begin{split}
							0\le x_{ij}^t \leq P_{ij} \cdot (Z_i-x_{ij}^{t-1})
						\end{split}
						\label{formula_4} %公式三引用标签， 通过 \eqref{formula_3}引用
					\end{equation}
				
				\item[(c)]三年种豆约束\par
				根据\eqref{formula_3}可确定当前可供耕地面积上界即最大值，现以参数 l 表示耕作时间（单位年）、以变量 t 表示耕作年份，再根据总耕地面积减去 l 年间所利用面积，可确定当前剩余最
				大土地面积即所取区间下界约束如下，后命名为三年种豆约束：
					\begin{equation}
						\begin{split}
							Z_i-\sum_{l=1}^{2} \sum_{j} x_{ij}^{t-l}\le x_{ij}^t\leq Z_i, ( j \in \textbf{豆类} , t \in (2024 , 2030] )
						\end{split}
						\label{formula_5} %公式四引用标签， 通过 \eqref{formula_4}引用
					\end{equation}
				\item[(d)] 水浇地种植面积约束\par
				水浇地种植时，水稻仅可一年种植一茬，而有其他作物可一年种植两茬，为便于计算，可将一年两茬作物视作拥有剩余耕地面积二倍的土地面积。故可确定
				水浇地种植面积约束如下：
					\begin{equation}
						\left\{\begin{array}{l}
							0\le 2 \cdot x_{i16}+\sum_{j=17}^{37}x_{ij}\leq 2 \cdot Z_{i} \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])   \\
							0\le\sum_{j = 17}^{34}x_{ij}\leq Z_{i}  \:, (0 \le x_{ij}\leq Z_i,i\in [27,34]) \\
							0\le\sum_{j = 35}^{37}x_{ij}\leq Z_{i}   \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])
						\end{array}\right.
						\label{formula_6} %公式五引用标签， 通过 \eqref{formula_5}引用
						\end{equation}
					\end{itemize}
			\end{samepage}

	根据上文，到此分别联立$\eqref{formula_1}$、$\eqref{formula_2}$、$\eqref{formula_3}$、$\eqref{formula_4}$、$\eqref{formula_5}$和$\eqref{formula_6}$可整理建立本题所需两大单目标优化数学模型如下：\par
	情况一：
	$\max \; Y_{t} =\sum_{i=1}^{54}\sum_{j = 1}^{41}(x_{ij}^{t} \cdot Q_{ij} \cdot pri_{ij}-x_{ij} \cdot C_{ij})-\sum_{j=1}^{41}max\left \{  \sum_{i=1}^{54}{x_{ij}^{t} \cdot Q_{ij}}-need_j,0\right \} \cdot pri_{ij}$
	$$
	s.t.
	\left\{\begin{array}{l}
		\sum_{j = 1}^{41} x_{ijk}^t\leq Z_{i} \: (k \in (1, 2), t \textgreater 2024) \\\\
		(Z_i-\sum_{l=1}^{2} \sum_{j} x_{ij}^{t-l})\le x_{ij}^t\leq Z_i, ( j \in \textbf{豆类} , t \in (2024 , 2030] ) \\\\
		0\le x_{ij}^t \leq P_{ij} \cdot (Z_i-x_{ij}^{t-1}) \\\\
		0\le 2 \cdot x_{i16}+\sum_{j=17}^{37}x_{ij}\leq 2 \cdot Z_{i} \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])   \\\\
		0\le\sum_{j = 17}^{34}x_{ij}\leq Z_{i}  \:, (0 \le x_{ij}\leq Z_i,i\in [27,34]) \\\\
		0\le\sum_{j = 35}^{37}x_{ij}\leq Z_{i}   \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])
	  \end{array}\right.
	$$
	  	\newline
	\begin{samepage}
		情况二：
		$\max \; Y_{t} =\sum_{i=1}^{54}\sum_{j = 1}^{41}(x_{ij}^{t} \cdot Q_{ij} \cdot pri_{ij}-x_{ij}^{t} \cdot C_{ij})-\sum_{j=1}^{41}max\left \{  \sum_{i=1}^{54}{x_{ij}^{t} \cdot Q_{ij}}-need_j,0\right \} \cdot pri_{ij} \cdot 50\%$
		$$
		s.t.
		\left\{\begin{array}{l}
			sum_{j = 1}^{41} x_{ijk}^t\leq Z_{i} \: (k \in (1, 2), t \textgreater 2024) \\\\
			(Z_i-\sum_{l=1}^{2} \sum_{j} x_{ij}^{t-l})\le x_{ij}^t\leq Z_i, ( j \in \textbf{豆类} , t \in (2024 , 2030] ) \\\\
			0\le 2 \cdot x_{i16}+\sum_{j=17}^{37}x_{ij}\leq 2 \cdot Z_{i} \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])   \\\\
			0\le x_{ij}^t \leq P_{ij} \cdot (Z_i-x_{ij}^{t-1}) \\\\
			0\le\sum_{j = 17}^{34}x_{ij}\leq Z_{i}  \:, (0 \le x_{ij}\leq Z_i,i\in [27,34]) \\\\
			0\le\sum_{j = 35}^{37}x_{ij}\leq Z_{i}   \:, (0 \le x_{ij}\leq Z_i,i\in [27,34])
		  \end{array}\right.
		$$
	\end{samepage} 

	\subsubsection{模型的求解}
	\begin{figure}[htbp]
		\centering
		\includegraphics[width=0.75 \textwidth]{photo_1.png} % 只需写文件名即可
		\caption{python代码表示图}
		\label{image_2} % 图片的引用标签， 通过 \ref{image_1}引用
	\end{figure}
	建立模型后，开始进行求解。本文主要采取python语言进行编程运算。将模型转化为python代码，以矩阵的形式进行运算求解得出最优结果如下（由于篇幅限制，仅展示2024年部分数据，其余结果见附件）：\par
		\begin{figure}[htbp]
			\centering
			\includegraphics[width=0.75 \textwidth]{excel_2.png} % 只需写文件名即可
			\caption{2024年度问题1超过部分浪费最优种植方案}
			\label{image_3} % 图片的引用标签， 通过 \ref{image_1}引用
		\end{figure}

		\begin{figure}[h]
			\centering
			\includegraphics[width=0.75 \textwidth]{excel_3.png} % 只需写文件名即可
			\caption{2024年度问题1超过部分降价出售最优种植方案}
			\label{image_4} % 图片的引用标签， 通过 \ref{image_1}引用
		\end{figure}

	\begin{samepage}
		\subsection{问题二模型的建立与求解}
		本题要求在结合考虑现实中可能发生市场变化不确定性和潜在种植风险的前提下，确定乡村2024$\sim$2030年农作物最优种植方案 
		\subsubsection{模型的建立：}
			根据题中所给信息可知，农作物种植成本平均每年增长5\%;
			不同作物售价变化率不同，粮食类作物售价稳定，蔬菜类年增长5\%， 食用菌年下降5\%

		综合上述内容，决定采用蒙特卡洛模拟的方法建立模型。
		由于希望得到最大化期望收益，同时考虑收益反差（其余风险指标）可将目标函数设定为：
		\begin{equation}
			\begin{split}
				\operatorname{Max} \; E(Z)-\lambda\cdot\operatorname{Var}(Z)
			\end{split}
			\label{formula_7}
		\end{equation}
		其中E(Z)为期望收益，$\lambda$为风险系数，Var(Z)为收益方差。
	\end{samepage}
	

	\subsubsection{模型的求解}
	建立模型后，开始进行求解。本文主要采取python语言进行编程运算。将题目一所构建模型结合蒙特卡洛模拟法进行运算求解得出最优结果如下（由于篇幅限制，仅展示2024年部分数据，其余结果见附件）：
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.75 \textwidth]{excel_4.png} % 只需写文件名即可
		\caption{2024年度问题2最优种植方案}
		\label{image_5} % 图片的引用标签， 通过 \ref{image_1}引用
	\end{figure}

	\begin{samepage}
		\subsection{问题三模型的建立与求解}
	问题三要求在问题二的基础上，进一步考虑作物之间的替代性和互补性，并且结合农作物预期销售量、销售价格和种植成本之间的相关性，为乡村制定出2024年到2030年的最优种植方案。
	\subsubsection{模型的建立}
	\textbf{1. 目标函数的确定}\par
	在问题三中，我们需要在问题二的基础上引入作物替代性和互补性。作物之间的替代性指的是当一种作物需求或价格下降时，可以用其他作物来替代；互补性指的是两种作物可以同时种植并带来收益提升。
	
	由此
	\begin{equation}
		\max Y = \sum_{i=1}^{n} \sum_{j=1}^{m} x_{ijk} \cdot (Q_{ij} \cdot pri_{ij} - C_{ij}) 
		+ \sum_{i=1}^{n} \sum_{j=1}^{m} \sum_{l=1, l \neq j}^{m} \left( \text{Comp}(x_{ijk}, x_{ikl}) - \text{Sub}(x_{ijk}, x_{ikl}) \right)
		\end{equation}
	其中，$\text{Comp}(x_{ijk}, x_{ikl})$ 和 $\text{Sub}(x_{ijk}, x_{ikl})$ 分别为互补性收益和替代性损失。

	\end{samepage}
	
	\textbf{2. 约束条件}\par
	\begin{itemize}
		\item[(a)] 土地面积约束\par
		每块土地上的种植面积不能超过该地块的可用面积：
		\begin{equation}
			\sum_{j=1}^{m} x_{ijk} \leq Z_i, \quad \forall i
		\end{equation}
		\item[(b)] 作物轮作约束\par
		每块地三年内必须至少种植一次豆类作物：
		\begin{equation}
			\sum_{t=1}^{3} x_{ijk} \geq M, \quad \forall j \in \text{豆类作物}
		\end{equation}
		\item[(c)] 替代性约束
		当作物 $j$ 和 $l$ 之间存在替代性时，种植作物 $j$ 的面积减少将导致作物 $l$ 的面积增加：
		\begin{equation}
			\text{Sub}(x_{ijk}, x_{ikl}) = \alpha \cdot (x_{ijk} - x_{ikl})
		\end{equation}
		\item[(d)] 互补性约束当作物 $j$ 和 $l$ 之间具有互补性时，两者种植面积增加将产生协同效应：
		\begin{equation}
			\text{Comp}(x_{ijk}, x_{ikl}) = \beta \cdot x_{ijk} \cdot x_{ikl}
		\end{equation}
	\end{itemize}
	

	\subsubsection{模型的求解}

	对于问题三的求解，我们将采用协方差矩阵和蒙特卡洛模拟的方法。具体步骤如下：
	\begin{enumerate}
		\item 协方差矩阵计算\par
		\hspace{2em}通过历史数据计算作物预期销售量、价格和种植成本的协方差矩阵，描述这些变量之间的相关性。
		\item 蒙特卡洛模拟\par
		\hspace{2em}使用蒙特卡洛模拟来模拟不同情况下的作物替代和互补效应，找到最优种植组合。
		\item 优化求解\par
		\hspace{2em}在考虑替代性和互补性的基础上，利用线性规划或其他优化方法求解最优种植方案。
	\end{enumerate}


	\section{模型的评价、改进与推广}
    \subsection{模型优点}
    \begin{enumerate}
        \item 优化模型灵活性\par
        \hspace{2em}本模型在处理不同作物的种植决策时，考虑了作物间的替代性与互补性，能够灵活应对农作物市场的需求变化。该模型允许在特定条件下，
        调整作物种植面积以最大化整体收益，并有效规避市场需求变化带来的风险。

        \item 数据驱动的决策支持\par
        \hspace{2em}通过引入销售价格、种植成本、产量等实际数据，本模型能够为决策者提供数据驱动的优化种植方案。这种方式能较为准确地反映现实情况，
        在不同的市场条件下，给出切实可行的种植建议。

        \item 考虑了多种约束条件\par
        \hspace{2em}模型不仅考虑了作物种植面积的约束，还引入了作物轮作、豆类作物的种植周期等多重约束条件，确保了种植方案的科学合理性。
    \end{enumerate}
    
    \subsection{模型缺点}
    \begin{enumerate}
        
        \item 替代性与互补性系数的估计\par
        \hspace{2em}模型中作物间的替代性与互补性系数采用了随机生成的方式，虽然在一定程度上能够模拟作物间的关系，但这些系数的估计可能缺乏实际依据。如果能基于更多的历史数据或试验结果来确定这些系数，模型的准确性将会大大提升。

        \item 对市场需求变化的敏感性\par
        \hspace{2em}模型对未来市场需求的预测基于当前的预期值，未充分考虑极端市场变化的情况。如果市场需求波动剧烈，模型可能无法快速响应并进行调整。
    \end{enumerate}
    
    \subsection{模型的改进}
    \begin{enumerate}
        
        \item 替代性与互补性系数的改进\par
        \hspace{2em}未来模型可以通过分析作物间的历史数据、市场反应，以及试验结果，更准确地估计替代性与互补性系数，从而提高模型的实际应用价值。这些系数可以通过机器学习等技术进行精确拟合，以提高模型的预测能力。

        \item 风险管理的引入\par
        \hspace{2em}在模型中引入更多的风险管理因素，例如利用蒙特卡洛模拟进一步评估市场需求和气候波动的风险，提供更加稳健的种植方案。同时，可以考虑添加更多关于农作物价格波动、气候变化等因素的情景分析，以增强模型对复杂现实的适应性。
    \end{enumerate}

	\newpage
	%Part Seven
	\section{参考文献}
	\vspace{-2em} % 减小上面的间距
	\begin{thebibliography}{9}  
		\bibitem{ref1} J. Kennedy and R.C. Eberhart. Particle swarm optimization. In Proceedings of the IEEE
		International Joint Conference on Neural Networks, pages 1942–1948, 1995.
		\bibitem{ref2} Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and et al. A fast and elitist multiobjective
		genetic algorithm: Nsga-ii. IEEE Transactions on Evolutionary Computation, 6(2), 2002.
		\bibitem{ref3} 司守全 孙玺菁，数学建模算法与引用，国防工业出版社 2021
		\bibitem{ref4} 姚志 and 何蒲明. 中国农村居民粮食消费需求及弹性测算. 统计与决策, 36(03):52–56,
		2020.
		\bibitem{ref5} 毛莉莎. 供应链视角下蔬菜批发市场定价策略及产销模式研究. 2023.
	\end{thebibliography}
	
	\newpage
	\section*{附录}
	
	附录1：支撑材料的文件列表
	\begin{table}[h]
    \center
	\begin{adjustbox}{max width=\textwidth} % 调整表格宽度为文本宽度
		\begin{tabular}{c c c}
			\hline
			文件名      & 含义            &      简介          \\
			\hline
			dataGandler.py       & python代码文件     &     清洗筛选数据      \\
			fun.py       & python代码文件     &      工具函数库      \\
			model.py       & python代码文件     &      数学模型    \\
			$is_single.csv$ & csv文件     &      初始化数据      \\
			Priceij.csv & csv文件     &      价格矩阵      \\
			Pij.csv & csv文件     &      可种植矩阵      \\
			Needj.csv & csv文件     &      需求矩阵      \\
			Cij.csv & csv文件     &      成本矩阵      \\
			Xij-2023.csv & csv文件     &      需求量矩阵      \\
			Zi.csv & csv文件     &      面积矩阵      \\
			Qij.csv & csv文件     &      产量矩阵      \\
			\hline
		\end{tabular}
	\end{adjustbox}
	\end{table}
 
	\newpage
	附录2：初始化代码和数据处理代码
	\begin{lstlisting}[language=python,columns=fullflexible,frame=shadowbox]
		import pandas as pd
		import fun

		type_to_symbol = {
			'平旱地': 'A',
			'梯田': 'B',
			'山坡地': 'C',
			'水浇地': 'D',
			'普通大棚': 'E',
			'智慧大棚': 'F'
		}
		symbol_to_type = {
			'A': '平旱地',
			'B': '梯田',
			'C': '山坡地',
			'D': '水浇地',
			'E': '普通大棚',
			'F': '智慧大棚'
		}

		# 新生成的表格
		new_table = pd.DataFrame()

		# 读取数据
		data1 = pd.read_csv('资料/2023年的农作物种植情况.csv', encoding='gbk')
		data1 = pd.DataFrame(data1[0:], columns=['种植地块', '作物编号', '作物名称', '作物类型', '种植面积', '种植季次'])

		data2 = pd.read_csv('资料/2023年统计的相关数据.csv', encoding='gbk')
		data2 = pd.DataFrame(data2[0:], columns=['序号', '作物编号', '作物名称', '地块类型',
												'种植季次', '亩产量', '种植成本', '销售单价'])
		new_table['作物名称'] = data2['作物名称']

		data3 = pd.read_csv('资料/乡村种植的农作物.csv', encoding='gbk')
		data3 = pd.DataFrame(data3[0:], columns=['作物编号', '作物名称', '作物类型', '种植耕地'])

		data4 = pd.read_csv('资料/农村的现有耕地.csv', encoding='gbk')
		data4 = pd.DataFrame(data4[0:], columns=['地块名称', '地块类型', '地块面积'])

		# 计算每亩收益
		shouyi = [fun.get_income(a, b, c) for a, b, c in zip(data2['亩产量'], data2['种植成本'], data2['销售单价'])]
		new_table['每亩收益'] = shouyi

		# 计算总产量
		zongliang = [0] * len(data2['作物名称'])

		for idx, name1 in enumerate(data2['作物名称']):
			for name2, s in zip(data1['作物名称'], data1['种植面积']):
				if name1 == name2:
					zongliang[idx] = zongliang[idx] + s
		new_table['总产量'] = zongliang

		arr = data3['种植耕地'].values
		for idx, e in enumerate(data3['种植耕地']):
			if type(e) is float:
				arr[idx] = data3['种植耕地'][idx - 1]
		data3['种植耕地'] = arr
		# print(data3)
		for name, di in zip(data3['作物名称'], data3['种植耕地']):
			col = [0] * len(data4['地块类型'])
			in_ = str(di).split('\r\n')
			for idx, e in enumerate(data4['地块类型']):
				if fun.in_di(in_, e):
					col[idx] = 1
			data4[name] = col

		data4.to_csv('./某个作物在某个地上是否能种植.csv')
		new_table.to_csv("./每亩收益和总产量.csv")

		# 获取Pij: 表示i号地中j号作物是否可以种植(1,0)
		Pij_org = pd.read_csv('某个作物在某个地上是否能种植.csv', encoding='utf-8', index_col=0)
		del Pij_org['地块类型']
		del Pij_org['地块面积']
		del Pij_org['地块名称']
		Pij_org.to_csv('数据/Pij.csv')
		# 作物列表
		P_rows = Pij_org.index
		P_cols = Pij_org.columns

		# Xij_2023: i号地j号作物种植情况
		Xij = pd.DataFrame(index=P_rows, columns=P_cols)


		Xij[:] = 0

		for idx, di_num, num, s in zip(range(len(data1['种植地块'])), data1['种植地块'], data1['作物编号'], data1['种植面积']):
			Xij.loc[data4['地块名称'].tolist().index(di_num), Xij.columns[num - 1]] = s
			# for a, b, c in zip(data2['作物编号'], data2['地块类型'], data2['亩产量']):
			#     if a == num and symbol_to_type[di_num[0]].strip() == b.strip():
			#         Qij.loc[data4['地块名称'].tolist().index(di_num), Qij.columns[num - 1]] = c
		Xij.to_csv('数据/Xij-2023.csv')

		# Qij: i号地j号作物产量
		Qij = pd.DataFrame(index=P_rows, columns=P_cols)
		Qij[:] = 0
		for num, ty, c in zip(data2['作物编号'], data2['地块类型'], data2['亩产量']):
			s = type_to_symbol[ty.strip()]
			for i in range(len(Qij)):
				if s in data4['地块名称'][i]:
					# print(i, data4['地块名称'][i], Qij.columns[num - 1], c)
					Qij.loc[i, Qij.columns[num - 1]] = c

		Qij.to_csv('数据/Qij.csv')

		# Priceij: i号地j号作物成本
		Priceij = pd.DataFrame(index=P_rows, columns=P_cols)
		Priceij[:] = 0
		for num, ty, sp in zip(data2['作物编号'], data2['地块类型'], data2['销售单价']):
			for idx in range(len(P_rows)):
				if type_to_symbol[ty.strip()] in data4['地块名称'][idx]:
					Priceij.loc[idx, Priceij.columns[num - 1]] = fun.get_singlePrice(sp)

		Priceij.to_csv("数据/Priceij.csv")

		# Cij: i号地j号作物成本
		Cij = pd.DataFrame(index=P_rows, columns=P_cols)
		Cij[:] = 0
		for num, ty, cost in zip(data2['作物编号'], data2['地块类型'], data2['种植成本']):
			for idx, name in zip(range(len(P_rows)), data4['地块名称']):
				if type_to_symbol[ty.strip()] in name:
					Cij.loc[idx, Priceij.columns[num - 1]] = cost

		Cij.to_csv('数据/Cij.csv')

		# Needj_2023: 作物j在2023年的需求量
		Needj = pd.DataFrame(columns=P_cols, index=['1'])
		for col in Qij.columns:
			Needj[col] = sum(Qij[col])
		Needj.to_csv('数据/Needj.csv')

		# Zi: 表示i号地的面积
		Zi = pd.DataFrame(columns=data4['地块名称'].tolist())
		Zi.loc[len(Zi)] = data4['地块面积'].tolist()
		Zi.to_csv('数据/Zi.csv')

		Zon = Qij * Xij * Priceij
		# print(Zon.sum().sum())

		is_single = pd.DataFrame(columns=['是否单季'])
		is_single['是否单季'] = [False] * len(data3['作物编号'])
		for idx, ji in zip(data2['作物编号'], data2['种植季次']):
			if '单季' == ji.strip():
				is_single.loc[idx - 1, '是否单季'] = True
		is_single.to_csv('数据/is_single.csv')
	\end{lstlisting}
	
	\begin{lstlisting}[language=python,columns=fullflexible,frame=shadowbox]
		import pandas as pd
		import numpy as np

		def get_income(a, b, c):
			arg1, arg2 = [float(arg) for arg in c.split('-')]
			return a * (arg1 + arg2) / 2 - b

		def get_singlePrice(c):
			arg1, arg2 = [float(arg) for arg in c.split('-')]
			return (arg1 + arg2) / 2


		def in_di(arr, arg):
			for e in arr:
				if arg in e:
					return True
			return False


		is_single = pd.read_csv('数据/is_single.csv', index_col=0)


		# 返回这块地占用的总面积（考虑单双季）
		def occ_area(z):
			s = 0
			for i, x in enumerate(z):
				s += x
				if not is_single.loc[i, '是否单季']:
					s -= x / 2
			return s

		def generate_matrix_combinations(Z, rows, cols):
			max_indices = np.ceil(Z / 0.1).astype(int)
			indices = np.zeros((rows, cols), dtype=int)
			while True:
				current_matrix = indices * 0.1
				yield current_matrix
				for i in range(rows):
					for j in range(cols):
						if indices[i][j] < max_indices[i] - 1:
							indices[i][j] += 1
							break
						indices[i][j] = 0
					else:
						continue
					break
				else:
					break
	\end{lstlisting}

	\begin{lstlisting}[language=python,columns=fullflexible,frame=shadowbox]
		import random

		import numpy as np
		import pandas as pd
		import copy
		import fun
		import itertools

		# 定义常量
		ZoneMax = 54
		CropMax = 41
		Is_Bean_Cons = False

		# 参数调整
		num_samples = 10  # 模拟次数
		loss_degree = 0  # 亏损程度enum(0, 0.5)
		sale_p = 0  # 小麦、玉米，预期销售量增长率 5% ~ 10%
		ge_sale_p = 0  # 其他作物预期销量增长率 -5% ~ 5%
		yield_P = 0  # 亩产量变化 -10% ~ 10%
		cost_p = 0  # 种植成本 5 %左右
		dish_price_p = 0  # 蔬菜类作物增长5%左右
		bact_price_p = 0  # 食用菌 1% ~ 5%
		t_price_p = 0.05  # 羊肚菌 5%
		lam = 0.5  # 风险厌恶系数

		# 读入数据
		X_t, P, C_t, Q_t, Z, Price_t, Need_t, X_Pre, X_PrePre = [0, 0, 0, 0, 0, 0, 0, 0, 0]


		def read_data():
			global X_t, P, C_t, Q_t, Z, Price_t, Need_t, X_Pre, X_PrePre
			X_t = pd.read_csv('数据/Xij-2023.csv', index_col=0)  # 表示i号地中j号作物在第t年的种植面积
			P = pd.read_csv('数据/Pij.csv', index_col=0)  # 表示i号地中是否可以种植j号作物
			C_t = pd.read_csv('数据/Cij.csv', index_col=0)  # 表示i号地中j号作物在第t年的种植成本
			Q_t = pd.read_csv('数据/Qij.csv', index_col=0)  # 表示i号地中j号作物在第t年的产量
			Z = pd.read_csv('数据/Zi.csv', index_col=0).iloc[0]  # 表示i号地的面积
			Price_t = pd.read_csv('数据/Priceij.csv', index_col=0)  # 表示i号地中j号作物在第t年的单价
			Need_t = pd.read_csv('数据/Needj.csv', index_col=0).iloc[0]  # 表示j号作物在t年的需求量
			X_Pre = pd.DataFrame(copy.copy(X_t))
			X_PrePre = pd.DataFrame(copy.copy(X_t))


		is_single = pd.read_csv('数据/is_single.csv', index_col=0)


		def sale_rate_change():
			Need_t.loc['小麦'] = (1 + sale_p) * Need_t.loc['小麦']
			Need_t.loc['玉米'] = (1 + sale_p) * Need_t.loc['玉米']
			for label in Need_t.index:
				if label in ['小麦', '玉米']:
					continue
				Need_t.loc[label] = (1 + ge_sale_p) * Need_t.loc[label]


		def yield_rate_change():
			for i in range(ZoneMax):
				Q_t.iloc[i].apply(lambda x: x * (1 + yield_P))


		def cost_rate_change():
			for i in range(ZoneMax):
				C_t.iloc[i].apply(lambda x: x * (1 + cost_p))


		def price_rate_change():
			for idx in range(CropMax):
				if 16 <= idx < 37:
					Price_t.iloc[:, idx].apply(lambda x: x * (1 + dish_price_p))
				if 37 <= idx < 40:
					Price_t.iloc[:, idx].apply(lambda x: x * (1 + bact_price_p))
				if idx == 40:
					Price_t.iloc[:, idx].apply(lambda x: x * (1 + t_price_p))


		def change_total():
			sale_rate_change()
			yield_rate_change()
			cost_rate_change()
			price_rate_change()


		# 返回这块地占用的总面积（考虑单双季）
		def occ_area(z):
			s = 0
			for i, x in enumerate(z):
				s += x
				if not is_single.loc[i, '是否单季']:
					s -= x / 2
			return s


		# 总面积约束
		def area_cons():
			a = [fun.occ_area(X_t.iloc[i]) for i in range(ZoneMax)]
			for idx, s, z in zip(range(len(a)), a, Z):
				if s > z:
					return False
			return True


		# 重茬种植约束
		def re_cons():
			for i in range(ZoneMax):
				for j in range(CropMax):
					up = P.iloc[i, j] * (Z.iloc[i] - X_Pre.iloc[i, j])
					if not (0 <= X_t.iloc[i, j] <= up):
						return False
			return True


		# 三年种豆约束
		def bean_cons():
			if not Is_Bean_Cons:
				return True
			for i in range(ZoneMax):
				down = Z.iloc[i] + sum(
					[X_Pre.iloc[i, j] + X_PrePre.iloc[i, j] for j in itertools.chain(range(0, 5), range(16, 19))])
				val = sum([X_t.iloc[i, j] for j in itertools.chain(range(0, 5), range(16, 19))])
				if not (down <= val <= Z.iloc[i]):
					return False
			return True


		# 水浇地约束1
		def water_area_cons1():
			for i in range(26, 34):
				val = X_t.iloc[i, 15] + sum([X_t.iloc[i, j] for j in range(16, 37)])
				up = 2 * Z.iloc[i]
				if not (0 <= val <= up):
					return False
			return True


		# 水浇地约束2
		def water_area_cons2():
			for i in range(26, 34):
				val = sum([X_t.iloc[i, j] for j in range(16, 34)])
				if not (0 <= val <= Z.iloc[i]):
					return False
			return True


		# 水浇地约束3
		def water_area_cons3():
			for i in range(26, 34):
				val = sum([X_t.iloc[i, j] for j in range(34, 37)])
				if not (0 <= val <= Z.iloc[i]):
					return False
			return True


		def cons_total():
			if area_cons() and re_cons() and bean_cons() and water_area_cons1() and water_area_cons2() and water_area_cons3():
				return True
			return False


		def Y():
			global loss_degree
			loss_degree = 0
			a = (X_t * Q_t * Price_t - X_t * C_t).sum().sum()
			b = 0
			for j in range(CropMax):
				for i in range(ZoneMax):
					b += max(X_t.iloc[i, j] * Q_t.iloc[i, j] - Need_t.iloc[j], 0)
			return a - b * (1 - loss_degree)


		income = []


		def model(num):
			read_data()
			income.clear()
			global X_t, Is_Bean_Cons, X_Pre, X_PrePre
			gen = fun.generate_matrix_combinations(Z.values, ZoneMax, CropMax)
			for i in range(2024, 2030):
				mx = 0.0
				if i >= 2025:
					Is_Bean_Cons = True
				for matrix in gen:
					temp = X_t
					X_t = pd.DataFrame(matrix)
					if not cons_total():
						X_t = temp
						continue
					mx = max(mx, Y())
				X_t.to_csv('./' + str(i) + '-' + str(num) + '.csv')
				X_PrePre = X_Pre
				X_Pre = X_t
				income.append(mx)


		# 最优方案索引
		g_idx = 0
		g_val = 0


		def mont_calo():
			global sale_p, ge_sale_p, yield_P, cost_p, dish_price_p, bact_price_p, t_price_p, g_idx
			for i in range(num_samples):
				sale_p = random.uniform(0.05, 0.1)
				ge_sale_p = random.uniform(-0.05, 0.05)
				yield_P = random.uniform(-0.1, 0.1)
				cost_p = random.uniform(0.045, 0.055)
				dish_price_p = random.uniform(0.045, 0.055)
				bact_price_p = random.uniform(0.01, 0.05)
				model(i)
				mean = np.mean(income)
				var = np.var(income)
				target = mean - lam * var
				if g_val > target:
					g_idx = i



	\end{lstlisting}
\end{document}